function [param,glob] = setup(param,glob,options)

    %% Create state space for customer base
    minb            = glob.minb;
    maxb            = glob.maxb;

 %   bgrid           = nodeunif(glob.n(1),(minb)^glob.curv,(maxb)^glob.curv).^(1/glob.curv);
 
    ubar = log(1 + log(1 + maxb-minb));
    ugrid = nodeunif(glob.n(1), 0, ubar);
    bgrid = minb + exp(exp(ugrid)-1)-1;


    bgrid0          = bgrid;

    
    %% Create idiosyncratic productivity matrix and grid
    [Pa, agrid] = setup_rouwen(param.rho_s, 0, param.sigma_s/(1 - param.rho_s^2), glob.n(2)-1);
    Pa = Pa';
    agrid = exp(agrid);  

    %% Find stationary distribution across productivities (non-superstar firms)
   
    a_stat = ones(size(agrid));
    a_stat = a_stat/sum(a_stat);
    a_stat = Pa'^10000*a_stat;


    %% Find stationary probabilities of being normal or superstar
   
    agrid   = [agrid;param.ebar];
    agrid   = log(agrid);
    agrid0  = agrid;

    % Update P
    Pa = [Pa*(1-param.p);a_stat'*(1-param.q)];
    Pa = [Pa,ones(glob.n(2),1)*param.p];
    Pa(end,end) = param.q;

    %% Find stationary distribution across productivities (all firms)
   
    a_stat = ones(size(agrid));
    a_stat = a_stat/sum(a_stat);
    a_stat = Pa'^10000*a_stat;


    %% Create finer grid for the histogram
    
    bgridf           = nodeunif(glob.nf(1),(minb)^glob.curv,(maxb)^glob.curv).^(1/glob.curv);
    Nbf              = size(bgridf, 1);

    sf               = gridmake(bgridf, agrid);
    Nsf              = size(sf, 1);
    
    glob.bgridf      = bgridf;
    glob.sf          = sf;
    glob.Nsf         = Nsf;
    
    glob.Bf = reshape(glob.sf(:,1), glob.nf(1), glob.nf(2));
    glob.Sf = reshape(glob.sf(:,2), glob.nf(1), glob.nf(2));

    
            
    %% Function space and nodes (fspace adds knot points for cubic splines)
    fspace          = fundef({'spli',bgrid,0,glob.spliorder(1)},...
                             {'spli',agrid,0,glob.spliorder(2)});

    sgrid           = funnode(fspace);
    s               = gridmake(sgrid);
    Ns              = size(s,1);

    %% Find actual grid given that fspace has added points for cubic spline

    %pgrid           = s(s(:,2)==s(1,2),1);
    %sgrid           = s(s(:,1)==s(1,1),2); % This is picking one of the other state and then 

    Nb              = size(bgrid,1);
    Ne              = size(agrid,1);

    glob.B = reshape(s(:,1), glob.n(1), glob.n(2));
    glob.S = reshape(s(:,2), glob.n(1), glob.n(2));
    
    %% Compute expectations matrix
    Phi             = funbas(fspace,s);
    PnK             = kron(Pa,speye(Nb));    % This is for expected value fn
    PnKf            = kron(Pa,speye(Nbf));
 
%    PnK             = kron(speye(Nb),Pa);    % This is for expected value fn
%    PnKf            = kron(speye(Nbf),Pa);


    %% Compute QZ matrix for approximating stationary distribution
    Qe              = kron(Pa,ones(Nb,1));
    
    Qef             = kron(Pa,ones(Nbf,1));

    glob.QeI             = kron(eye(size(Pa)),ones(Nbf,1));

    % construct Q matrix for iterating on productivity but not labor
    fspaceerg   = fundef({'spli',glob.bgridf,0,1});    
    Qb          = funbas(fspaceerg, sf(:,1));

    glob.QeNl     = dprod(Qef,Qb);

    %% 
     fspaceerg   = fundef({'spli',agrid,0,1});
     entrant     = agrid - param.d_E;
     entrant     = min(entrant, max(agrid));
     entrant     = max(entrant, min(agrid));
     Qshift      = funbas(fspaceerg, entrant);

    glob.a_stat_E = a_stat;
    glob.a_stat_E(end) = 0; % ensure no entrants are superstars upon entry
    glob.a_stat_E = glob.a_stat_E/sum(glob.a_stat_E);
    
    glob.a_stat_E = Qshift'*glob.a_stat_E;
    
    glob.a_stat_E = glob.a_stat_E';     


    %% Transition matrix from signal to realized productivity
     
    % just get a matrix that redistributes mass from the finer signal grid
    % to the coarser productivity grid
    
%    fspaceerg   = fundef({'spli',agrid,0,1});
%    entrant     = sig_grid;
%    glob.Qshift      = funbas(fspaceerg, entrant');
    
%    glob.sig_grid = sig_grid';
    
    %% evenly distribute mass across labor choices
     G               = zeros(Nsf,1);

    for i = 1:glob.n(2)
        idx = find(sf(:,2) == agrid(i));
        
        G(idx) = glob.a_stat_E(i)/length(idx);
    end

    glob.entrant_prod_dist = G;
    
    %% get distribution of entrants over both states now
        % Distribution of entrant firms
           
    %% Create other one time only basis matrices
    Phi2            = splibas(agrid0,0,glob.spliorder(2),s(:,2));
%    Phi3            = splibas(zgrid0,0,glob.spliorder(3),s(:,3));

%%
    glob.ve_l = zeros(length(glob.a_stat_E), 1);


    %% Declare additional global variables
    glob.bgrid0     = bgrid0;
    glob.bgrid      = bgrid;

    glob.sgrid0     = agrid0;
    glob.sgrid      = agrid;

%     glob.zgrid0     = zgrid0;
%     glob.zgrid      = zgrid;

    glob.Ps          = Pa;
%     glob.Pz          = Pz;

    glob.Qe         = Qe;

    glob.Ne         = Ne;
    glob.Np         = Nb;

    glob.fspace     = fspace;
    glob.s          = s;
    glob.Ns         = Ns;
    glob.Nb         = Nb;
    
    glob.Phi2       = Phi2;
%    glob.Phi3       = Phi3;

    glob.Phi        = Phi;
    glob.PnK        = PnK;
    glob.PnKf       = PnKf;
    
    glob.agrid = agrid;
    glob.a_stat = a_stat;
    
    glob.Qef        = Qef;
    
end
